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Abstract 


Epigenetic variation allows for rapid changes in phenotypes without 
alterations to nucleotide sequences. These epigenetic signatures may di- 
verge over time among isolated populations. Epigenetic incompatibility 
following secondary contact between these populations could result in the 
evolution of reproductive isolating mechanisms. If epigenetic incompati- 
bility drove the evolution of species isolating mechanisms, we expect to see 
significant epigenetic differentiation between these species. Alternatively, 
epigenetic variation could be the result of predominantly environmental 
variables and not align along species boundaries. A methylation sensitive 
amplified fragment length polymorphism analysis was performed on in- 
dividuals of the closely related katydid species Neoconocephalus robustus 
and N. bivocatus. We observed significant variation in total methylation 
levels between species. However, genetic differentiation remained larger 
than epigenetic differentiation between species groups. We measured a sig- 
nificant correlation between the epigenetic and genetic distance between 
individuals. Epigenetic differentiation is therefore likely the result of an 
interaction between genetic and epigenetic loci and not a mechanism for 
species differentiation. We therefore did not find evidence to support our 
hypothesis of an epigenetically mediated mechanism for speciation be- 
tween N. robustus and N. bivocatus. 
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Introduction 


Epigenetic variation can lead to changes in phenotypic expres- 
sion without any change to the nucleic acid sequence (Berger et 
al. 2009). Unlike genetic variation, epigenetic changes can pro- 
duce reversible, heritable phenotypic changes within a lineage 
Jablonka and Lamb 1998). Epigenetically controlled phenotypic 
plasticity may therefore play an important role in rapid, adaptive 
changes (Led6on-Rettig et al. 2012, Burggren and Crews 2014). 
DNA methylation, the most commonly studied form of epige- 
netic modification, involves the addition of methyl groups, usu- 
ally to CpG dinucleotides that regulate gene expression (Boyes 
and Bird 1991). Methylation regulated phenotypic plasticity has 


been documented across multiple taxonomic groups (Braam and 
Davis 1990, Rapp and Wendel 2005, Varriale and Bernardi 2006, 
Elango et al. 2009, Herrera and Bazaga 2010). Until recently, epi- 
genetically dependent variation has been largely overlooked as an 
evolutionary mechanism contributing to speciation (Smith et al. 
2016). Facultative phenotypes can be maintained in a population 
as adaptive alternatives to divergent environmental conditions. 
These alternative phenotypes may ultimately become fixed in a 
population by way of genetic assimilation (Pal and Miklos 1999, 
West-Eberhard 2005). 

Phenotypic variation can also evolve in response to epigenetic 
incompatibility (= hybrid incompatibility caused by epialleles; 
Bente and Scheid 2017, Blevins et al. 2017) between groups. Fol- 
lowing the epigenetic diversification of groups in isolation, inter- 
mediate forms may show reduced fitness following secondary con- 
tact Jablonka and Lamb 1999). Aphids moved to a new host plant 
and allowed to reproduce asexually quickly developed a preference 
for the new host and produced morphological characteristics simi- 
lar to conspecifics utilizing this host. Backcrossing to the parental 
line resulted in the production of non-viable offspring (Shaposh- 
nikov 1966). Heritable phenotypic variation in aphids has been 
correlated with changes in methylation (Field et al. 1989, Walsh 
et al. 2010). Methylation dependent incompatibility could be rein- 
forced by reproductive isolation (Pal and Miklds 1999). 

In many animals, behavioral isolation plays a significant role 
in maintaining species boundaries. Acoustic communication has 
been studied as a mechanism for reproductive isolation in both 
vertebrate and invertebrate groups (Gerhardt and Huber 2002). 
Neoconocephalus, a New World genus of katydids, possesses a di- 
verse range of habitat preferences (Walker et al. 1973, Frederick 
2013) and call phenotypes (Schul and Patterson 2003, Beckers and 
Schul 2008, Deily and Schul 2009). Females utilize these calls in 
mate recognition and in phonotaxis, the directional movement to- 
wards a calling male (Greenfield 1990, Triblehorn and Schul 2009, 
Beckers and Schul 2010). The acoustic communication system of 
Neoconocephalus allows for reproductive isolation among the mul- 
tiple species that may be found in sympatry (Schul et al. 2014). 

The sibling species (Snyder et al. 2009) N. robustus (Scudder, 
1862) and N. bivocatus Walker, Whitesell, & Alexander, 1973 dif- 
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fer both in their habitat preferences and call types. Neoconocepha- 
lus robustus is a grassland generalist and utilizes grasslands with a 
wide range of flora and environmental conditions. Neoconoceph- 
alus bivocatus, on the other hand, is a prairie obligate preferring 
drier habitats composed primarily of tall prairie grasses (Walker 
et al. 1973). Both species co-occur in some areas, e.g. at the edges 
of prairie remnants or along creeks within prairies. This species 
pair is largely morphologically cryptic but maintains divergent 
call characters and call preferences (Deily and Schul 2004, 2006). 
While the two species can be distinguished using genetic mark- 
ers (e.g. mitochondrial haplotypes or AFLP markers (Snyder et 
al. 2009, Ney and Schul 2017)), genetic differentiation is weak 
and evidence for hybridization (e.g. intermediate calls) occur fre- 
quently. Population-level genetic variation was shown to be very 
low, with almost no geographically associated genetic variation 
between sites as distant as 448 km apart (Ney and Schul 2017). 
This suggests that these two species diverged very recently, most 
likely during the current glacial cycle (Frederick 2013). 

A likely scenario of this divergence event includes an ancestral 
population living in both habitat ranges. In each habitat, different 
epigenetic patterns would be expressed, which then became fixed 
within each subpopulation, resulting in decreased hybrid fitness. 
Genetic differences leading to the differences in male calls and fe- 
male preferences evolved later, in part driven by reproductive rein- 
forcement (Deily and Schul 2004). Thus we hypothesize here that 
epigenetic differentiation, driven by the different environmental 
conditions experienced in their core habitats, drove divergence of 
the species. The lack of substantial genetic differentiation between 
N. robustus and N. bivocatus could point to an underlying epige- 
netic incompatibility driving species differentiation. This hypoth- 
esis predicts that the epigenetic differentiation between the species 
would be larger and more stereotyped than the genetic differences. 

DNA methylation has been described in many insect groups 
including multiple orthopteran species (Sarkar et al. 1992, Robin- 
son et al. 2011). Methylation sensitive amplified fragment length 
polymorphism (MS-AFLP) analysis is one of the few techniques 
that allows for the quantification of genome-wide patterns of cy- 
tosine methylation without any previous knowledge of genome 
sequences (Xiong et al. 1999). The MS-AFLP technique uses 
two isoschizomer restriction enzymes (MspI and Hpall) that 
recognize the same restriction site, cutting in the same location 
(5’-CACGG-3’), but have different sensitivities to the presence 
of cytosine methylation. By comparing the presence/absence of 
fragments produced by both enzymes, the methylation status at 
each restriction site can be evaluated. In addition, genetic poly- 
morphisms can be evaluated between individuals that lack the 
restriction site in both digestions as an indication of a change in 
nucleotide sequence. 

To distinguish whether divergence of these two species was 
initiated by epigenetic, rather than genetic variation, we addressed 
two questions using a MS-AFLP technique. First, we asked to what 
extent epigenetic and/or genetic differentiation predicted species 
assignment. Second, we analyzed whether the patterns of genome- 
wide DNA methylation were correlated with genetic variation. 


Materials and methods 


Specimen collection.—We utilized a total of 94 males collected in 
the summers of 2006, 2013, and 2014 from eight grassland sites 
around the state of Missouri (Suppl. material 1: Table S1). We 
used males’ calls to localize individuals in the field and collect- 
ed them by hand after sunset. We identified males in the field as 
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belonging to members of the target group through their call and 
morphological features, including cone pigmentation and body 
size, prior to collection (Walker et al. 1973). As the two species 
are morphologically cryptic, we identified individuals as either N. 
bivocatus or N. robustus based on their divergent call characteristics 
as described below. Individuals were collected from among eight 
grassland sites within Missouri and from within three different 
years (Suppl. material 1: Table S1). Thus, a range of epigenetic 
patterns, expressed within variable environmental conditions, 
was included in each species sample. Fixed differences observed 
between the two species are therefore likely to be species specific 
rather than the result of differences in environmentally mediated 
methylation patterns. 


Call recordings. —We recorded male calls in 2006 and 2013 within 
three days of collection using an Audiotechnica ATR 55 micro- 
phone and a Marantz PMD-671 solid-state recorder (16 bit, 48 
kHz sampling rate). Recordings were made outdoors with males 
placed in individual mesh cages (approximately 10x20x10 cm) 
spaced at least 3 m apart. In 2014, we recorded male calls in the 
field immediately preceding collection using a Tascam DR-40 lin- 
ear PCM recorder (16 bit, 48 kHz sampling rate). Ambient tem- 
peratures ranged from 22-28°C. 


Temporal call analysis. —Call temporal patterns were analyzed and 
species assignments made as described in Ney and Schul (2017). 
In short, we marked each sound pulse produced during a single 
closing movement of the wings (Walker 1975) using custom soft- 
ware. The recordings were rectified and the envelope extracted 
with a temporal resolution of 0.125 ms. Pulses were detected 
automatically using a threshold-based algorithm and manually 
checked before saving the data to text-files. We analyzed about 2 
s of each male’s call containing 150-250 pulses. Further analyses 
were conducted on the text files in MS Excel. 

Male calls of N. robustus have a single pulse rate of about 200/s, 
equivalent to a pulse period of about 5 ms (Walker 1975). Females 
recognize this pattern by the absence of silent gaps longer than 
about 2 ms (Deily and Schul 2004). In contrast, male calls of N. 
bivocatus have two alternating pulse periods of about 4.5 ms and 
7 ms (Deily and Schul 2004), resulting in a ‘galloping rhythm’ 
or ‘pulse pairs’. Female N. bivocatus recognize this pattern by the 
rate of the pulse pairs: pulse pair rates around 87/s are attractive 
and attractiveness decreases towards faster and slower rates. For 
the pulse pairs to be detectable, the two alternating pulse periods 
must differ sufficiently (J. Schul, unpublished work). 

We quantified the ratio of the means of the alternating pulse 
periods (longer pp / shorter pp). For the single pulse calls of N. ro- 
bustus, this should result in values close to one, while in N. bivoca- 
tus significantly larger values result (Bush and Schul 2010). Among 
the calls of 94 males analyzed, there were two distinct groups of 
values, one between 1 and 1.18 and a second group with ratios 
>1.38. A natural break occurred in the data between 1.18 and 
1.38. We classified ratios <1.18 as ‘N. robustus’, and >1.38 as ‘N. 
bivocatus’ in call type. 


Spectral call analysis. —Call spectra were analyzed and species as- 
signments made as described in Ney and Schul (2017). In short, 
we analyzed the amplitude spectrum of male calls using a fast Fou- 
rier transformation (FFT, Hamming window, frame length 256) in 
Audacity v.1.3 (Audacity Team 2008). Spectra were averaged over 
2 s of the call. The main energy in the spectrum of Neoconocephalus 
calls is concentrated in a narrow band between 7 and 15 kHz with 
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the peak frequency differing among species (Schul and Patterson 
2003). We measured the center frequency of this low frequency 
band as the geometric mean of the upper and lower cut-off fre- 
quencies at -3 dB from the peak amplitude. 

In N. bivocatus, the center frequencies of this band have a mean 
of about 10 kHz among individuals, ranging from 7 to 15 kHz. Fe- 
males have little spectral selectivity in this frequency range (Deily 
and Schul 2006). The center frequencies of calls with a pulse peri- 
od ratio of >1.38 that also fell into this range of center frequencies 
were Classified as N. bivocatus. One individual had a pulse period 
ratio greater than 1.38 and a center frequency lower than 7 kHz; 
this male was classified as an intermediate caller, as its call charac- 
teristics differed from that of both species. 

In N. robustus, the low frequency band is narrower than in N. 
bivocatus and is typically limited to 10 kHz and below (Schul and 
Patterson 2003). Indeed, frequencies well above 10 kHz have an 
inhibitory effect on female phonotaxis (Deily and Schul 2006). Of 
the individuals with pulse period ratios <1.18, center frequencies 
of all but five individuals clustered between 6 and 9 kHz. The five 
remaining individuals had center frequencies ranging from 9.7 to 
11.1 kHz, which would be less attractive for N. robustus females. 
These individuals were classified as possessing an intermediate 
call type. All intermediate callers were removed from the analysis 
of epigenetic and genetic differentiation. 


Molecular analysis.—We aim here to analyze the differences of 
trans-generational methylation patterns between two species. 
Therefore, we selected tissue that is likely to have little tissue-spe- 
cific methylation differences between the two species. We collected 
tissue from hind femurs (mostly muscle and cuticle) that should 
be under similar selection in both species and thus little differen- 
tial tissue-specific methylation. Using leg tissue also avoids the risk 
of contamination from GI tract microbiota. 

We removed the hind femurs of collected males and placed 
them in 95% EtOH for DNA preservation. We later extracted DNA 
from the hind femurs using the DNeasy Blood & Tissue Kit (Qia- 
gen Inc., Valencia, CA, USA). DNA quantification was performed 
on each sample by spectrophotometry (NanoDrop 1000, Thermo 
Scientific, Wilmington, DE). Genomic DNA was stored at -80°C 
prior to molecular analysis. 

We used a MS-AFLP assay modified from Xu et al. (2000). 
DNA (55 ng) from each sample was digested in two separate dou- 
ble digest reactions (EcoRI/HpalI and EcoRI/MsplI). EcoRI selec- 
tively cuts the sequence 5’-GAATTC-3”. Hpall and MspI are isoschi- 
zomers, meaning they selectively cut at the sequence (5’-CCGG- 
3’), but differ in their sensitivity to cytosine methylation at those 
sites. Hpall will only cut if the external cytosine is hemimethyl- 
ated. MspI cuts when the internal cytosine is either hemi or fully 
methylated. Both enzymes will cut if the target site is completely 
unmethylated. Using this method, we evaluated the CpG methyla- 
tion of restriction sites by comparing the fragments produced by 
both digests. 

Digestion and ligation were carried out together to prevent 
regeneration of restriction sites. Synthetic double stranded DNA 
adaptors (Xu et al. 2000) were ligated to the cleaved ends of re- 
striction sites. The EcoRI/Hpall and EcoRI/MspI digestion/ligation 
reaction (11 ul final volume) is comprised of 1.1 pL 10X CutSmart- 
™ buffer (New England Biolabs, USA), 0.55 pg/pL of Bovine Se- 
rum Albumin (BSA) (New England Biolabs), 0.3 pL water, 5 U of 
EcoRI HE, 1 U of either Hpall or MspI, 1 U T4 ligase (New England 
Biolabs), 1 pL ATP (10 pM), EcoRI adaptors (5 pM), either Hpall 
or MspI adaptors (50 pM), and 5.5 pL genomic DNA (10 ng/ pL). 
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Table 1. Loci produced by selective primer combinations used in 
the MS-AFLP analysis. Shown are the primer pair combinations, 
the number of scored bands per primer pair, and the number of 
those bands classified as polymorphic methylation sensitive loci 
(MSL), and polymorphic genetic loci. 


EcoRI Msp1/Hpall Bands MSL Genetic loci 
-AAC -ATC 277 183 265 
-AGC -AAT 22 49 102 
-AGC -ATC 318 170 301 

Total 822 402 668 


The reaction was incubated at 37°C for 2 hours. Preselective PCR 
was conducted with 1:10 dilute digestion/ligation products, the 
Eco+A (5’-GACTGCGTACCAATTCA-3’), and the HpallI/MspI+A 
(5’-GATGAGTCTIAGAACGGA-3’) primers using thermocycler set- 
tings as described in Snyder et al. (2009). We performed selective 
PCR independently with three primer pairs (Table 1) and 1:100 
dilute preselective PCR products from both the Hpall and MspI 
reactions. Fluorescently labeled Eco primers (Eco+AAC (6GFAM), 
Eco+AGC (PET)) were used in selective PCR (as described in Sny- 
der et al. 2009) and the products multiplexed and diluted to pro- 
duce a 1:10 final dilution of each product. Fragments were separat- 
ed in an ABI 3730 genetic analyzer at the DNA Core Facility, Uni- 
versity of Missouri. We called MS-AFLP bands using GeneMarker 
v.1.6 (Hulce et al. 2011) using an automated peak-calling scheme 
(as described in Holland et al. 2008) and called alleles between 
75-500 bp with a minimum peak intensity of 50. 


Data analysis. —We obtained presence/absence fragment data for 
both EcoRI/Hpall and EcoRI/MspI datasets from GeneMarker 
(Hulce et al. 2011). Presence of a fragment in both the EcoRI/ 
Hpall and EcoRI/MspI digestions (1/1) was described as unmeth- 
ylated. The presence of a fragment in one digestion but not the 
other (0/1 or 1/0) was defined as methylated (either hemimethyl- 
ated or internal cytosine methylated, respectively). If fragments 
were absent in both digestions (0/0) the loci were considered un- 
informative as this state could be due to either full methylation at 
the target site (methylation of both the inner and outer cytosine) 
or the absence of the fragment due to variation in the nucleotide 
sequence between individuals (Schulz et al. 2013, Fulneéek and 
Kovarik 2014). 

The methylation sensitivity of each locus was identified us- 
ing the MSAP Package (Pérez-Figueroa 2013) in R v.3.1.2 (R Core 
Team 2015). Loci were classified as methylation sensitive loci 
(MSL) if there was evidence of methylation in at least 5% of the 
sampled individuals at that locus. Genetic markers were extracted 
from the MS-AFLP loci. Fragments that were present in one or 
both MspI and Hpall analyses were scored as present. Fragments 
that were absent in both reactions were scored as absent. This 
method of using MS-AFLP loci to estimate genetic parameters has 
been found to produce similar results to standard AFLP markers 
(Smith et al. 2016). Only polymorphic MSL and genetic loci that 
differed among sampled individuals were used in further analy- 
ses. We evaluated variation for each epigenetic (MSL) and genetic 
locus individually then calculated the mean diversity of MSL and 
genetic loci using Shannon’s diversity index (I). We compared the 
relative frequency of total methylation (internal cytosine meth- 
ylated and hemimethylated fragments) and unmethylated MSL 
between species. Significant variation between species was tested 
using a Mann-Whitney U test. 
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We tested for differences in MSL and genetic differentiation us- 
ing two-way analyses of molecular variance (AMOVA; Excoffier et 
al. 1992) that grouped individuals again by species. Significance 
of the test statistic’s (O,,) deviation from zero was estimated based 
on 1000 permutations of individuals among groups. Principal co- 
ordinate analyses (PCoA) of both epigenetic and genetic loci were 
performed using the R stats package v.3.1.2 (cmdscale; R Core 
Team 2015) to visualize the Euclidean distance between species. A 
population structure analysis was implemented for all individuals 
using the program STRUCTURE v.2.3.3 (Pritchard et al. 2010). The 
admixture model was used, allele frequencies correlated, with a 
run length of 100,000 (Burnin = 10,000) for 10 replicates each of K 
= 1-10 (genetic clusters). The estimate of the most well supported 
K was calculated as described in Evanno et al. (2005) and imple- 
mented in Structure Harvester v.0.6.94 (Earl 2012). The program 
Clumpp v.1.1.2 (Jakobsson and Rosenberg 2007) was used to align 
the 10 repetitions of the most well supported number of clusters. 

If similar signals exist for the MSL and genetic structure then 
epigenetic and genetic distance may be significantly correlated. 
We estimated the Euclidean distance between individuals for both 
epigenetic and genetic datasets using the R stats package v.3.1.2 
(R Core Team 2015) and compared the distance matrices using a 
Mantel test (Mantel 1967; ade4 package v.1.7.2; Dray et al. 2007). 
We utilized 10,000 permutations of the design matrix to deter- 
mine the significance of the correlation coefficient. 


Results 


Genome-wide variation in methylation.—We included 94 males in 
our analysis. We determined the males’ species assignments based 
on the temporal and spectral frequency call preferences of N. ro- 
bustus and N. bivocatus. The call analysis classified the 94 males as 
31 N. bivocatus, 57 N. robustus, and 6 having an intermediate call 
type; five of the intermediates had the N. robustus temporal pat- 
tern and one the N. bivocatus pattern, with frequency spectra that 
fell outside of these respective species’ specific pattern (Fig. 1). In- 
dividuals possessing intermediate call phenotypes were removed 
from further analyses. The MS-AFLP analysis yielded 227, 277, and 
318 fragments in each of three selective PCR reactions (Table 1). Of 
these fragments, a total of 364 were polymorphic for their meth- 
ylation status among sampled individuals and were classified as 
MSL. A loci-calling scheme, utilizing the MS-AFLP dataset, allowed 
for the inference of the genetic state of fragments among individu- 
als, similar to a traditional AFLP analysis. In total, 668 polymor- 
phic genetic loci were produced from 822 MS-AFLP fragments. 

We compared whether species differed significantly in their 
proportion of genome-wide methylated sites. The relative frequen- 
cy of genome-wide methylation (combined hemimethylation and 
internal cytosine methylation) showed low (Fig. 2; N. robustus = 
0.272; N. bivocatus = 0.240) but significant variation between spe- 
cies (Mann-Whitney U, W = 587, p = 0.001243). 


Epigenetic and genetic structure.—We examined epigenetic and ge- 
netic diversity as they relate to species assignment. The within-spe- 
cies epigenetic Shannon diversity index was 5.0616 + 0.2054 and 
5.0094 + 0.2087, within N. robustus and N. bivocatus, respectively. 
N. robustus and N. bivocatus’ genetic Shannon diversity indexes 
were 5.2802 + 0.2242 and 5.2179 + 0.1943, respectively. Diver- 
sity in the genetic loci, as measured in this study, was significantly 
greater than that of the MSL (epigenetic) diversity measured in 
both species (Wilcoxon rank sum test; N. robustus, W = 2760, p < 
0.0001; N. bivocatus, W = 779, p = 0.0001). 
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Fig. 1. Species assignment based on call pulse period ratio and 
center frequency. Labeled boxes indicate the calls classified as N. 
robustus and N. bivocatus. Individuals that fall outside of species 
classifications were removed from further epigenetic and genetic 
analyses (as described in Ney and Schul 2017). 
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Fig. 2. Comparison of genome-wide methylation levels between 
species. Mann-Whitney U test; p<0.005 (**). Between-species sig- 
nificant variation is in total methylation (internal cytosine methyl- 
ated and hemimethylated fragments). 


Epigenetic and genetic structure was evaluated using a two- 
level AMOVA. The between-species epigenetic divergence was ®,, = 
0.0504 (P < 0.0001). This is slightly less than one-third of the genet- 
ic divergence observed between species, O,, = 0.1591 (P < 0.0001). 
Greater genetic variation between species would suggest that genetic 
mechanisms are underlying species differentiation (Table 2). 
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Table 2. Two-level AMOVA of MSL or genetic loci produced from 
MS-AFLP markers among populations grouped by species as- 
signment. Included are genetic and epigenetic variance between 


groups, ®,,, and the corresponding p value. 
Between group Within group ©,  pvalue 
variance variance 
MSL (epigenetic) 3.429 (5.04%)  64.54(94.96%) 0.0504 <0.0001 


Genetic loci 15.38 (15.91%) 81.28 (84.09%) 0.1591 <0.0001 


The principal coordinate analyses (PCoA), calculated using MSL 
and genetic profiles, showed that genetic variance was smaller than 
epigenetic variance within species (Fig. 3). The epigenetic PCoA 


A 


C1 (11.4%) 


15 


between species elucidated little meaningful differentiation among 
groups, with ellipses (95% confidence intervals) overlapping al- 
most entirely. The PCoA of genetic data between species explained 
17.7% of the genetic variance in the first two coordinates (Fig. 3A) 
and showed moderate variation between individuals based on call 
phenotype, with species showing divergence in Euclidean space. 
Three N. robustus individuals clearly fell into the N. bivocatus cluster 
within the genetic PCoA (Fig. 3A), indicating a possible mismatch 
between phenotypic and genotypic assignments. 

The Bayesian analysis of genetic structure revealed the best- 
supported number of genetic clusters to be K = 2, based on AK 
values (Fig. 4B; Evanno et al. 2005). These two genetic clusters 
aligned closely to species boundaries (Fig. 4A). As in the PCoA, 
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Fig. 3. PCoA of N. robustus and N. bivocatus utilizing genetic (A) and epigenetic (B) data. Plotted are the two most informative principal 
components calculated for the genetic and epigenetic loci datasets, as derived from the MS-AFLP fragment analysis. A. Genetic Euclid- 
ean distance with individuals grouped by species assignment. B. Epigenetic Euclidean distance with individuals grouped by species as- 
signment. Group labels show the centroid of the points for each group. The long axis of the ellipse represents the direction of maximum 


dispersion and the short axis the direction of minimum dispersion. 
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Fig. 4. Consensus shared ancestry population structure for epigenetic and genetic loci. A. and C. Bar plots using MS-AFLP loci to estimate ge- 
netic (A) and epigenetic (C) structure among N. robustus and N. bivocatus using the software package STRUCTURE. B. and D. Delta K graphs 
for K = 1-10 genetic clusters showing moderate support for K = 2 genetic clusters (B) and low support for K = 4 epigenetic clusters (D). 
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Fig. 5. Scatterplot of between-individual Euclidean genetic and 
epigenetic distance showing significant positive correlation be- 
tween genetic and epigenetic differentiation. The correlation was 
tested using a Mantel test and 10,000 permutations of the design 
matrix to determine significance. 


the same three individuals identified as N. robustus were assigned 
primarily to the N. bivocatus genetic cluster, indicating that they 
possess the N. robustus call type but a N. bivocatus genotype. As in 
the PCoA, the epigenetic structure analysis did not show support 
for significant epigenetic (MSL) structure between species groups. 
AK values did not rise above 30 for any permutation in the num- 
ber of clusters (Fig. 4D), an indication of low support for epige- 
netic structure. In addition, the best supported number of clusters, 
K = 4, showed no grouping of epigenetic clusters by species assign- 
ment (Fig. 4C). 

While significant epigenetic variation was observed between 
species, significant genetic variation was also detected. We inves- 
tigated the correlation between inter-individual genetic and epi- 
genetic Euclidean distance. Between-individual epigenetic and ge- 
netic distance showed a strong positive correlation (Fig. 5, Mantel 
R = 0.8448, p = 0.0001). 


Discussion 


We found significant differentiation among species in both 
epigenetic and genetic markers. Genetic differentiation, however, 
was larger than epigenetic differentiation between species. The 
AMOVA and PCoA analyses indicated that methylation patterns 
varied among individuals but showed little differentiation be- 
tween species. Epigenetic distance among individuals correlated 
with inter-individual genetic differentiation, suggesting that the 
low levels of epigenetic differentiation between species have been 
pulled along by genetic differentiation. 


Genetic differentiation.—Genetic differentiation between species 
was low but significant, as would be expected between two closely 
related taxa (Snyder et al. 2009, Ney and Schul 2017). Genetic 
structure aligned with species boundaries, with the exception of 
the three phenotypically N. robustus individuals that showed a 
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majority assignment to the primarily N. bivocatus genotypic clus- 
ter. These mismatch genotype/phenotype assignments were con- 
firmed for these three individuals in both the genetic structure 
analysis (Fig. 4A) as well as the genetic PCoA (Fig. 3A). Similar 
occurrences of mismatched genotype/phenotype individuals have 
been found between N. robustus and N. bivocatus previously (Fred- 
erick 2013, Ney and Schul 2017). These mismatched genotype/ 
phenotype individuals are likely the result of recent hybridiza- 
tion events. A severe drought during the sampling period (Ney 
and Schul 2017) led locally to sharp population declines of one 
or both species. Both reductions in population sizes and environ- 
mental disturbances can increase rates of hybridization between 
closely related species (Lamont et al. 2003, Seehausen 2006). 


Epigenetic variation.—N. robustus showed a significantly higher lev- 
el of gnome wide methylation than N. bivocatus (Fig. 2). This may 
be due to variation in the species’ habitat preference. Variable en- 
vironmental conditions can cause shifts in epigenetics, observed 
between diverging natural habitats as well as between natural and 
altered habitat types (Gao et al. 2010). N. robustus, unlike N. bivo- 
catus, is a grassland generalist and may show epigenetic patterns 
resulting from exposure to a more variable set of environmental 
conditions during a lifetime and across generations. Because of 
this, N. robustus may possess a larger repertoire of adaptive genes 
for the varied environmental conditions it may utilize. The higher 
level of methylation found in N. robustus may therefore be the re- 
sult of higher levels of methylation required to silence the large 
repertoire of adaptive genes when not in use. 

While epigenetic differentiation between species was signifi- 
cant, it remained lower than genetic differentiation. This study did 
not show support for the epigenetic regulation of species-specific 
phenotypes; however, this does not eliminate a possible epige- 
netic mechanism underlying phenotypic differentiation. While 
the MS-AFLP technique has many benefits, there are also some 
inherent limitations to its application. For example, MS-AFLPs can 
underestimate genome-wide levels of methylation (Fulneéek and 
Kovarik 2014). As a genome-wide scanning method MS-AFLP only 
detects methylation at 5'-CCGG-3' sites. In addition, it is unable to 
discriminate between full methylation states (hypermethylation 
of both cytosines) and sequence variation at or near the restric- 
tion site (Xiong et al. 1999). Phenotypic variation controlled by 
a smaller number of methylated sites is unlikely to be detected. A 
similar investigation into gnome-wide methylation found no sig- 
nificant differentiation in methylation between the solitarious and 
gregarious phases of Locusta migratoria (Robinson et al. 2011), de- 
spite divergent phenotypes and the identification of differentially 
expressed genes between morphs (Kang et al. 2004). Thus, this 
technique is most valuable as a first step in investigating genome- 
wide methylation patterns. 


Correlation between epigenetic and genetic diversity.—MS-AFLP vari- 
ation often shows correlations with genetic variation (Liu et al. 
2012). In our study, individuals showing greater genetic variation 
on average showed greater epigenetic variation as well (Fig. 5). 
Changes in DNA methylation over time are correlated with ge- 
netic relatedness and suggest that DNA methylation maintenance 
may be under genetic control (Bjornsson et al. 2008). In addition, 
genetic variation in retrotransposons (i.e. their presence/absence) 
could affect the methylation state of retroelements, requiring more 
methylation if more repetitive elements are present (Michaud et 
al. 1994). Repetitive elements have been estimated to make up 
thirty percent of the genome of the orthopteran L. migratoria (Wil- 
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more and Brown 1975). In addition, the retrotransposon SINE is 
differentially methylated between the solitarious and gregarious 
phases of the species (Guo et al. 2010). 

Mechanisms of neutral evolution could also account for the 
correlation between epigenetic and genetic variation. In the event 
of substantial gene flow between species, strong divergent selec- 
tion would be needed to maintain divergent epigenetic variation 
between species. Gene flow between groups would reduce differ- 
entiation accumulated via drift. Evidence from this study, however, 
suggests that genetic differentiation between species is significant 
and gene flow therefore relatively low (Table 2). Stochastic process- 
es of drift then could allow epigenetic patterns to diverge between 
species in parallel, resulting in correlations between epigenetic and 
genetic variation without a causal link (Richards et al. 2010). 


Conclusions 


We hypothesized that the phenotypic variation observed be- 
tween N. robustus and N. bivocatus was the result of an epigenetic- 
mediated mechanism of species differentiation that led to genetic 
isolation. While we found clear evidence of genomic methylation 
in Neoconocephalus, our findings did not support a role for meth- 
ylation in species isolation. Both the lower level of epigenetic dif- 
ferentiations and the correlation between inter-individual epige- 
netic and genetic diversity support the alternative hypothesis, i.e. 
that differences in methylation patterns between species evolved 
in response to genetic variation. Epigenetics may still play a key 
role in phenotypic differentiation within Neoconocephalus katydids 
through the differential regulation of a relatively small number of 
genes of large effect; a mechanism not detectable with the meth- 
ods used here. Further work identifying differentially expressed 
genes between call types could allow for the targeted analysis of 
methylation patterns at these sites. 
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